!***************************************************************************
MODULE DATA_ARRAY       
        real,dimension(2500,2500) :: Rij_n,Rij
        real,dimension(2500,2500) :: pairtype
        integer :: Qn,pairn
        integer,dimension(2500) :: pairatomi,pairatomj,pairresi,pairresj
        integer :: atomnum,resnum,matnum,ndom,ntype
END MODULE DATA_ARRAY
!****************************************************************************
! writed by wyong 2010.8.9
! updated by Lingci 2015/6/28

        program qnatcon

        use DATA_ARRAY

        implicit none

        integer :: i,j,k,atomi,atomj,resi,resj,ios,MODELI,coli,colj,cold
	integer :: Q,Qt,nargs,itype,pairnum
	real :: Qfrac,Qplu,Qhelix1,Qhelix2,ll
!	integer,dimension(20) :: connum,Qs,Qsfrac,Qsplu
	integer :: connum,Qs,Qsfrac,Qsplu
        character(len=99) :: line,listfile,datfile,arg
	character(len=8) :: ctmp,QMODEL
	character(len=20),dimension(20) :: col
	real :: dist
        real,dimension(8) :: tmp
        real,dimension(500,3) :: Coor
        character(len=50) :: HEADER
        character(len=3) :: atomtype

        nargs = iargc()
        if(nargs < 8)then
                write(*,*) "QCAL.x 1          2    3    4    5        6       7      8"
                write(*,*) "QCAL.x maplist    coli colj cold pdb_traj Atomnum Qmodel ll"
                write(*,*) "QCAL.x binding.list 1    2    3    tmp.pdb  196     Q      0.42"
                write(*,*) "Q algorithm: Q, Qfrac, Qplu"
                goto 1000
        elseif(nargs == 8)then
                call getarg(1,arg)
                listfile = arg

                call getarg(2,arg)
                read(arg,'(I2)') coli
                call getarg(3,arg)
                read(arg,'(I2)') colj
                call getarg(4,arg)
                read(arg,'(I2)') cold

                call getarg(5,arg)
                datfile = arg

                call getarg(6,arg)
                read(arg,'(I4)') atomnum
!                write(*,'(A20,I5)') 'Total atom number ',atomnum

                call getarg(7,arg)
                QMODEL=arg

                call getarg(8,arg)
                read(arg,'(F7.3)') ll
        endif

        !================================
	connum=0

        open(10,file=listfile,status='old')
        readlist: do
                read(10,('(A99)'),iostat=ios) line
                if(ios .ne. 0) exit readlist
		!read(line,'(I4,1X,I4,F7.3,1X,F5.2,1X,I2)') atomi,atomj,dist,ratei,itype
		connum=connum+1
		read(line,*) (col(i),i=1,5)
		read(col(4),'(I4)') resi
		read(col(5),'(I4)') resj
		read(col(coli),'(I4)') atomi
		read(col(colj),'(I4)') atomj
		read(col(cold),'(F7.3)') dist
!		write(*,'(I5,I5,I5,F7.3)') connum,atomi,atomj,dist
		Rij_n(atomi,atomj)=dist
		pairatomi(connum)=atomi
		pairatomj(connum)=atomj
		pairresi(connum)=resi
		pairresj(connum)=resj
	enddo readlist
        close(10)

	pairn=connum
!        write(*,'(A11,I5,A20)') 'There is ',pairn,' LJ pairs'

        !==================================
        ! open trajectory file in PDB format
        open(20,file=datfile,status='old')
	if ( QMODEL == 'Q' ) then
        open(30,file='Q.dat',status='unknown')
        open(40,file='Qnorm.dat',status='unknown')
	endif
        if ( QMODEL == 'F' ) then
        open(50,file='Qfrac.dat',status='unknown')
	endif
        if ( QMODEL == 'P' ) then
        open(60,file='Qplu.dat',status='unknown')
	endif
        if ( QMODEL == 'H' ) then
        open(70,file='Qhelix.dat',status='unknown')
        endif
        readtrj: do
                read(20,('(A70)'),iostat=ios) line 
                if(ios .ne. 0) exit readtrj
                read(line,*) HEADER

                if ( HEADER == 'MODEL' ) then
                        read(line,'(6X,I8)') MODELI
			if ( mod(MODELI,10000) == 0) then
			write(*,'(I6)') MODELI
			endif
                endif

                if ( HEADER == 'ATOM' ) then
                   read(line,'(13X,A2)') atomtype
                   if ( atomtype == 'CA' .or. atomtype == 'CP' .or. atomtype == 'CB' ) then
                        read(line,'(6X,I5)') atomi
                        read(line,'(30X,3F8.3)') (tmp(i),i=1,3)
                        Coor(atomi,1)=tmp(1)
                        Coor(atomi,2)=tmp(2)
                        Coor(atomi,3)=tmp(3)
                   endif
                endif

                if ( HEADER == 'ENDMDL' ) then
                        do i=1,atomnum
                           do j=1+1,atomnum
                                Rij(i,j)=sum( (Coor(i,1:3)-Coor(j,1:3))**2 )
                                Rij(i,j)=sqrt(Rij(i,j))
                           enddo
                        enddo

        		if ( QMODEL == 'Q' ) then
                        call qcalc(Q)
                        write(30,'(I7,1X,I7)') MODELI,Q
                        write(40,'(I7,1X,F7.3)') MODELI,Q*1.0/pairn
			endif
		        if ( QMODEL == 'F' ) then
                        call qcalc2(Qfrac)
                        write(50,'(I7,1X,F7.3)') MODELI,Qfrac*1.0/pairn
			endif
		        if ( QMODEL == 'P' ) then
                        call qcalc3(Qplu,ll)
                        write(60,'(I7,1X,F7.3)') MODELI,Qplu*1.0/pairn
			endif
                        if ( QMODEL == 'H' ) then
                        call qcalc4(Qhelix1,Qhelix2)
                        write(70,'(I7,1X,F7.3,1X,F7.3)') MODELI,Qhelix1*1.0/pairn,Qhelix2*1.0/pairn
                        endif
                endif

        enddo readtrj

        close(10)
        close(20)

        if ( QMODEL == 'Q' ) then
        close(30)
        close(40)
	endif
        if ( QMODEL == 'F' ) then
        close(50)
	endif
        if ( QMODEL == 'P' ) then
        close(60)
	endif
        if ( QMODEL == 'H' ) then
        close(70)
        endif

1000    end program

!========================================
        subroutine qcalc(Q)

        use DATA_ARRAY

        implicit none
        integer, intent(out):: Q 
        integer :: ipair,i,j,npair
        real :: R1,R0
	real,parameter :: RATE=1.2
        
        Q=0
	npair=pairn

        do ipair=1,npair
		i=pairatomi(ipair)
		j=pairatomj(ipair)
                R1=Rij(i,j)
                R0=Rij_n(i,j)
		if ( R1 <= (RATE * R0) ) then
                	Q=Q+1
                endif
        enddo

        end subroutine 

!===============================================
! Reference: Robert B. Best, Jeetain Mittal
! Balance Between alpha and beta structures in Ab initio Protein Folding
! JPCB, 2010, 114, 8790-8798
! Qaa=sum(1 + 1/exp(beta*(Rij-RATE*R0ij)) )
! beta= 5.0 (1/angstrom)
! lambda= 1.2
        subroutine qcalc2(Qfrac)

        use DATA_ARRAY

        implicit none
        real, intent(out):: Qfrac
        integer :: ipair,i,j,npair
        real :: R1,R0
	real,parameter :: Beta=5.0
	real,parameter :: RATE=1.2

        Qfrac=0.0
        npair=pairn

        do ipair=1,npair
                i=pairatomi(ipair)
                j=pairatomj(ipair)
                R1=Rij(i,j)
                R0=Rij_n(i,j)
		Qfrac = Qfrac + 1.0 / ( 1 + exp(Beta*(R1-RATE*R0)) )
        enddo

        end subroutine

!===============================================
! Reference: PLUMED: A portable plugin for free-energy calculations
! Sij=1; if R=Rij-D0<=0
! Sij= 1-(R/R0)^n / 1-(R/R0)^m ; if R=Rij-D0>0, m>n
! Default value: m=12, n=6, D0=2.0, R0=Rij_nat

        subroutine qcalc3(Qplu,l)

        use DATA_ARRAY
        
        implicit none
        real, intent(in) :: l
        real, intent(out):: Qplu
        integer :: ipair,i,j,npair
        real :: R1,R0,D0
        real,parameter :: m=20, n=10
                
        Qplu=0.0
        npair=pairn
            
        do ipair=1,npair
                i=pairatomi(ipair)
                j=pairatomj(ipair)
                R1=Rij(i,j)
                R0=Rij_n(i,j)
                D0=l*R0
                Qplu = Qplu + (1-((R1-D0)/R0)**n) / (1-((R1-D0)/R0)**m)
        enddo

        end subroutine

!========================================
!========================================
!========================================

        subroutine qcalc4(Qhelix1,Qhelix2)

        use DATA_ARRAY

        implicit none
        real, intent(out):: Qhelix1,Qhelix2
        integer :: ipair,i,j,resi,resj,npair,Helix(2000)
	real :: Qpair(2000)
        real :: R1,R0
        real,parameter :: Beta=5.0
        real,parameter :: RATE=1.2

        Qhelix1=0.0
        Qhelix2=0.0
        npair=pairn
	Qpair(:)=0.0

        do ipair=1,npair
                i=pairatomi(ipair)
                j=pairatomj(ipair)
                resi=pairresi(ipair)
                resj=pairresj(ipair)
                R1=Rij(i,j)
                R0=Rij_n(i,j)
                if ( R1 <= (R0 + 1.0) ) then
			Qhelix1=Qhelix1+1.0 / ( 1 + exp(Beta*(R1-RATE*R0)) )
!			Qpair(resi)=1
			Qpair(resi)=1.0 / ( 1 + exp(Beta*(R1-RATE*R0)) )
                endif
        enddo

        i=1
        do while ( i < 1000 )
		if ( Qpair(i)>0 .and. Qpair(i+1)>0 .and. Qpair(i+2)>0 ) then
			Helix(i)=1
			Helix(i+1)=1
			Helix(i+2)=1
		endif
		i=i+1
	enddo

	i=1
	do while ( i < 1000 )
		if (Helix(i)==1) then
			Qhelix2=Qhelix2+Qpair(i)
		endif
		i=i+1
	enddo

        end subroutine  

